import numpy as np
from sympy import *
import matplotlib.pyplot as plt

# 根据mathematica解出模型表达式
'''
    Solve[{
        Sin[\[Theta]] == 0.64*Sin[\[Beta]] + 0.77*Cos[\[Beta]] Cos[t], 
        \[Beta] == 23.5 Degree*Sin[(216/180)*Pi], 
        x0 == 3/Tan[\[Theta]]}, 
        {\[Beta], \[Theta], t, x0}]
'''
# 计算影长的函数
fx=lambda t:3*cot(asin(0.7477*cos(t)-0.1543))

# 画出曲线图
time=np.linspace(9,15,100)
x=np.array([float(abs(t-12)*15*pi/180) for t in time])
y=np.array([fx(xi) for xi in x])
plt.plot(time,y,color="red")
plt.title("Sun shadow length curve from 9:00 to 15:00 in Tiananmen Square")
plt.xlabel("time")
plt.ylabel("shadow length/(m)")
plt.savefig('problem_1.png')
plt.show()
